﻿using System;
using System.Collections.Generic;
using System.Text;
using BarStructures.Deformed;
using SimpleLinearAlgebra;
using MathNet.Numerics.LinearAlgebra;

namespace BarStructures.Analysers
{
    public class FDCondensation_XY_DM: Analyser
    {
        private Matrix eigenVectors = null;

        private Vector eigenValues = null;

        private double _ev;

        public FDCondensation_XY_DM(Model model, string name, double ev)
            : base(model, name)
        {
            this._ev = ev;
        }

        public override void Analyse()
        {
            DeformedState ds = new FrameDeformedState();
            ds.AssignToModel(model);

            if (!UseExistingIndexation)
                model.SetAutoSystemIndexes();

            int modelDimension = model.GetNodesCount() * 2;

            //подготовка матриц
            SimpleMatrix stiffness = new SimpleMatrix(modelDimension, modelDimension);
            SimpleMatrix mass = new SimpleMatrix(modelDimension, modelDimension);

            //построение матриц жесткости и масс
            ISimpleMatrix elementMatrix;
            int barsCount, ind1, ind2;
            barsCount = model.GetBarsCount();
            for (int i = 0; i < barsCount; i++)
            {
                Bar bar = model.GetBar(i);
                ind1 = bar.NodeOne.SystemIndex * 2;
                ind2 = bar.NodeTwo.SystemIndex * 2;

                elementMatrix = ElementXYHelper.BarGlobalStiffnessMatrix(bar);
                stiffness.AddSubMatrix(ind1, ind1, elementMatrix.GetSubMatrix(0, 0, 2, 2));
                stiffness.AddSubMatrix(ind1, ind2, elementMatrix.GetSubMatrix(0, 2, 2, 2));
                stiffness.AddSubMatrix(ind2, ind1, elementMatrix.GetSubMatrix(2, 0, 2, 2));
                stiffness.AddSubMatrix(ind2, ind2, elementMatrix.GetSubMatrix(2, 2, 2, 2));

                elementMatrix = ElementXYHelper.BarDescreeteMassMatrix(bar);
                mass.AddSubMatrix(ind1, ind1, elementMatrix.GetSubMatrix(0, 0, 2, 2));
                mass.AddSubMatrix(ind1, ind2, elementMatrix.GetSubMatrix(0, 2, 2, 2));
                mass.AddSubMatrix(ind2, ind1, elementMatrix.GetSubMatrix(2, 0, 2, 2));
                mass.AddSubMatrix(ind2, ind2, elementMatrix.GetSubMatrix(2, 2, 2, 2));
            }

            //DbgPrint.SimpleMatrix(mass, "mass_0.txt");
            //DbgPrint.SimpleMatrix(stiffness, "stiffness_0.txt");

            //индексы основных степеней свободы
            List<int> primaryDegreesList = new List<int>();
            int[] primaryDegrees;
            for (int i = 0; i < model.GetNodesCount(); i++)
            {
                Node node = model.GetNodeBySystemIndex(i);
                if (node.DegreeSelection.X)
                    primaryDegreesList.Add(i * 2);
                if (node.DegreeSelection.Y)
                    primaryDegreesList.Add(i * 2 + 1);
            }
            primaryDegrees = primaryDegreesList.ToArray();

            //индексы второстепенных степеней свободы
            List<int> secondaryDegreesList = new List<int>();
            int[] secondaryDegrees;
            for (int i = 0; i < model.GetNodesCount(); i++)
            {
                Node node = model.GetNodeBySystemIndex(i);
                if (!node.Support.Sx && !node.DegreeSelection.X)
                    secondaryDegreesList.Add(i * 2);
                if (!node.Support.Sy && !node.DegreeSelection.Y)
                    secondaryDegreesList.Add(i * 2 + 1);
            }
            secondaryDegrees = secondaryDegreesList.ToArray();

            //индексы для матрицы составленной из основных и второстепенных степеней свободы
            int[] mergedDegrees = new int[secondaryDegrees.Length + primaryDegrees.Length];
            Array.Copy(primaryDegrees, mergedDegrees, primaryDegrees.Length);
            Array.Copy(secondaryDegrees, 0, mergedDegrees, primaryDegrees.Length, secondaryDegrees.Length);

            //Получение подматриц
            ISimpleMatrix Krr = stiffness.GetSubMatrix(primaryDegrees, primaryDegrees);
            ISimpleMatrix Ksr = stiffness.GetSubMatrix(secondaryDegrees, primaryDegrees);
            ISimpleMatrix Krs = stiffness.GetSubMatrix(primaryDegrees, secondaryDegrees);
            ISimpleMatrix Kss = stiffness.GetSubMatrix(secondaryDegrees, secondaryDegrees);
            ISimpleMatrix Mrr = mass.GetSubMatrix(primaryDegrees, primaryDegrees);
            ISimpleMatrix Mss = mass.GetSubMatrix(secondaryDegrees, secondaryDegrees);
            ISimpleMatrix Mrs = mass.GetSubMatrix(primaryDegrees, secondaryDegrees);
            ISimpleMatrix Msr = mass.GetSubMatrix(secondaryDegrees, primaryDegrees);
            ISimpleMatrix Kmrg = stiffness.GetSubMatrix(mergedDegrees, mergedDegrees);
            ISimpleMatrix Mmrg = mass.GetSubMatrix(mergedDegrees, mergedDegrees);


            //-ev * Is
            ISimpleMatrix Is = new SimpleMatrix(secondaryDegrees.Length, -this._ev);
            //(Krr * Mrs + Krs * Mss) 
            ISimpleMatrix save1 = MxOp.Add(MxOp.Product(Krr, Mrs), MxOp.Product(Krs, Mss));
            //inv(Ksr * Mrs + Kss * Mss - ev * Is)
            ISimpleMatrix save2 = MxOp.Add(MxOp.Product(Ksr, Mrs), MxOp.Product(Kss, Mss));
            save2 = MxOp.Inverse(MxOp.Add(save2, Is));
            //(Ksr * Mrr + Kss * Msr)
            ISimpleMatrix save3 = MxOp.Add(MxOp.Product(Ksr, Mrr), MxOp.Product(Kss, Msr));
            //всё вместе
            ISimpleMatrix D = MxOp.Product(MxOp.Product(MxOp.Product(save1, save2), save3), -1);

            //Матрица  преобразования C = Krr * Mrr + Krs * Msr + D;
            ISimpleMatrix C = MxOp.Add(MxOp.Add(MxOp.Product(Krr, Mrr), MxOp.Product(Krs, Msr)), D);


            //решение конденсированной проблемы собственных значений
            Matrix ir_C = new Matrix(((SimpleMatrix)C).GetInternalData());

            EigenvalueDecomposition evd = new EigenvalueDecomposition(ir_C);
            eigenVectors = evd.EigenVectors;
            eigenValues = evd.RealEigenvalues;
        }

        public override DeformedState GetDeformedStateForEigenValue(double eigenValue)
        {
            return null;
        }

        public override List<double> GetEigenValuesList()
        {
            if (eigenValues == null)
                Analyse();
            List<double> result = new List<double>();
            for (int i = 0; i < eigenValues.Length; i++)
                result.Add(eigenValues[i]);
            return result;
        }

        public override List<double> GetEigenVecsListForEigenValue(double eigenValue)
        {
            return null;
        }

        public override List<double> GetFullEigenVecsListForEigenValue(double eigenValue)
        {
            throw new NotImplementedException();
        }
    }
}
